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V  ABSTRACT 


The  concept  of  a  locking  material  and  the  potential  effectiveness  of 
the  locking  material  for  possible  use  as  a  countermeasure  to  specified 
mine  blasts  have  been  investigated  in  this  report.  The  effectiveness  of 
the  countermeasure  has  been  demonstrated  by  the  use  of  analytical  methods 
and  second  order  accurate  computer  codes.  Computer  codes  have  been  devel¬ 
oped  to  investigate  and  evaluate  such  a  design  that  uses  the  counter¬ 
measures.  These  computer  codes  include  programs  in  one  space  dimension, 
and  axisymmetric  coordinates. 

As  a  result  of  the  investigation,  specific  design  parameters,  candi¬ 
date  materials  and  fields  for  further  study  have  been  identified. 
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\  SUMMARY  • 

The  concept  of  a  locking  material  and  the  potential  effectiveness  of  the  locking 
material  for  possible  use  as  a  countermeasure  to  specified  mine  blasts  have  been 
investigated  in  this  report.  The  effectiveness  of  the  countermeasure  has  been 
demonstrated  by  the  use  of  analytical  methods  and  second  order  accurate  computer 
codes.  Computer  codes  have  been  developed  to  investigate  and  evaluate  such  a  design 
that  uses  the  countermeasures.  These  computer  codes  include  programs  in  one  space 
dimension,  and  axisymmetric  coordinates. 

As  a  result  of  the  investigation,  specific  design  parameters,  candidate  materials 
and  fields  for  further  study  have  been  identified. 
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1.  INTRODUCTION 

This  report  presents  the  results,  to  date,  of  investigations  concerning  the 

development  of  a  countermeasure  to  mine  blasts.  In  particular,  the  investigations  have 

been  conducted  to  develop  computer  simulation  techniques  that  are  capable  of  simulating 
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the  transient  dynamic  response  of  Locking  materials  *  or  locking  material-structure 
combinations.  These  deveoioped  techniques  have  been  used  to 

(a)  examine  if  a  locking  material  shield  has  the  potential  of  protecting  a  given 
structure  from  dynamic  loads  that  have  a  time  history  similar  to  that  of  mine 
blasts,  and 

(b)  provide  a  tool  for  designing  the  locking  material  shield-structure 
combination. 
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Terms  such  as  foams,  distended  materials  or  nonreactive  porous  solids  ,  have  been 
used  to  describe  the  locking  materials.  Typically,  a  locking  material  is  characterized  by 
a  hydrostatic  pressure-density  curve  similar  to  that  shown  in  Figure  1.1.  The  model 
described  in  Figure  1.1  is  an  idealized  behavior.  In  practice,  a  more  realistic  model  may 
be  necessary.  A  material,  which  follows  the  model  shown  in  Figure  1.2  behaves  like  an 
elastic  solid  below  a  pressure  Pg.  For  pressure  p>  pg,  the  pore  spaces  collapse  and  the 
material  locks  at  a  density  .  The  subsequent  behavior  of  the  material  is  that  of  an 
incompressible  material.  In  practical  shock  interaction  and  attenuation  calculations  by 
numerical  techniques,  the  model  shown  in  Figure  1.1  is  very  often  replaced  by  a  pressure 
density  relationship  similar  to  those  shown  in  Figures  1.2  or  1.3.  In  Figure  1.3,  after 
reaching  the  locking  density,  the  material  is  assumed  to  follow  a  pressure-density 
behavior  of  the  corresponding  solid.  Furthermore,  unloading  will  follow  different  paths 
as  shown  in  the  Figures  1.2  and  1.3.  Such  unloading  paths  are  necessary  to  reflect  the 
fact  that  a  collapsed  porous  space  can  not  be  recovered.  An  exception  to  this  rule  is  the 
behavior  of  graphite  foam^. 

The  theoretical  foundation  for  expecting  a  locking  material  to  be  an  effective  peak 
stress  attenuator  or  a  countermeasure  to  mine  blasts  can  be  summarized  as  follows. 
Materials,  that  approximate  a  pressure-density  behavior  of  Figures  1.1,  1.2  or  1.3,  imply 
a  large  ratio  of  unloading  or  rarefaction  wave  velocity  to  the  loading  or  shock  velocity. 
The  large  ratio  is  the  result  of  the  fact  that  the  velocity  of  the  loading  shockwave  is 
related  to  the  slope  of  a  straight  line  connecting  the  initial  and  final  states  in  the 
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pressure  density  diagram.  The  velocity  of  the  unloading  wave  or  rarefaction  wave  is 
related  to  the  slope  of  the  pressure-density  curve  for  the  solid  at  the  locked  state.  The 
attenuation  is  caused  by  the  unloading  or  rarefaction  wave  overtaking  the  loading  shock 
wave  because  of  the  relatively  higher  velocity  of  the  unloading  waves.  In  practical 
situations,  the  effects  of  shear  stresses  and  yielding  need  to  be  considered. 

1.1  Candidate  Materials 


In  principle,  locking  material  can  be  produced  from  any  parent  solid  material.  The 
needed  process  of  production  requires  the  formation  of  a  nonreactive  porous  solid  of  a 
density  lower  than  that  of  the  parent  material.  A  commerical  aluminum  locking  material 
MO-AL  (Emerson  and  Cummings)  has  been  available  in  the  past.  This  is  usually 
inhomogeneous.  Homogeneous  aluminum  locking  material  can  be  produced  by 
techniques^  such  as  hot  pressing  of  aluminum  powders,  cold  pressing  of  aluminum 
powders  followed  by  sintering  or  a  repeated  sequence  of  hot  pressing  that  is  followed  by 
sintering  of  aluminum  powders  and  the  use  of  silica  microballons.  Commerically, 
graphite  locking  materials  are  available  from  companies  such  as  National  Carbon 
Company.  Of  course,  a  very  common  locking  material  of  comparatively  low  threshold 
pressure  capability  that  is  easily  available  is  the  styro-foam.  Depending  on  the 
environment  and  the  peak  stresses  that  need  to  be  attenuated,  different  types  of  locking 
materials  can  be  chosen  or  designed  for  a  particular  application.  For  example,  in  the 
case  of  a  tank  that  is  subjected  to  mine  blasts,  a  locking  material  made  of  steel  may  be 
an  interesting  possibility.  In  some  cases,  a  combination  of  different  locking  materials  or 
a  sandwich  construction  may  be  more  practical. 


1.2  Theoretical  Foundation 


First,  a  locking  material  in  a  state  of  one  dimensional  strain  along  x-axis  is 
considered.'  A  pressure-density  relationship  as  shown  in  Figure  1.4  has  been  assumed. 
Initially,  stresses  and  velocities  in  the  region  x>0,  have  been  assumed  to  be  equal  to  zero. 
At  time  t  =  0,  it  is  assumed  that  a  stress  -ox  =  a0  «  applied  at  the  left  boundary.  The 
left  boundary  is  initially  at  x  =  xQ  (t  =  0)  =  0  (figure  1.5).  In  addition  to  the  pressure 
density  behavior,  the  behavior  in  shear  is  assumed  to  be  elastic-plastic  with  the  following 
yield  condition. 
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(tJj-tJj)2  +  (c2-a3)2  +  (a3-Oj)2  =  2Y2  (1.1) 

In  this  equation,  o  j,  a2*  o3  are  the  principal  stresses  and  Y  is  the  yield  stress  in 
simple  tension.  For  the  assumed  one-dimensional  strain  conditions,  a  a  x  -p  relationship 
can  be  derived  from  the  pressure-density  relationship  and  the  assumed  yield  condition. 
The  resulting  ox  -p  relationship  is  illustrated  in  Figure  1.6.  The  slope  of  the<jx  -p  curve 
for  p  2  p^is  derived  from  the  mechanical  behavior  of  the  solid  that  constitutes  the  parent 
material  for  the  given  locking  material  of  initial  density  pQ.  For  example,  p  for  an 
aluminum  locking  material  can  be  in  the  range  of  1.1  gms/c.c.  to  2.3  gms/c.c.  The  value 
of  is  approximately  equal  to  2.7  gms/c.c.  In  practice,  p^  is  usually  slightly  below  the 
value  of  the  solid  density  of  the  parent  material. 

For  t>0,  the  applied  stress  at  the  left  boundary  results  in  two  stress  waves 
propagating  in  the  positive  x-direction.  The  first  wave  is  an  elastic  forerunner  that 
carries  the  stress  discontinuity  that  corresponds  to  Y*.  The  second  wave  is  a  shock  wave 
that  carries  the  stress  discontinuity  that  corresponds  to  the  difference  between  oq  and 
Y*.  The  shockwave  velocity  is  denoted  by  k .  The  following  jump  conditions  are  valid 
across  the  elastic  forerunner  and  the  shock  wave. 


P„<ke  -  V  -  P,«<.  -  V 

(1.2) 

A*  •  • 

Y  =  -  pfi(ke  -  vft)  ve  +  PQ(ke  -  vq)  vq 

(1.3) 

P*(W  *  pe(kr  ve} 

(1.4) 

a  -  Y  =  -  v^)  v4  +  p  e(k4  -  vg)  ve 

(1.5) 

In  these  equations,  v  and  v  are  the  particle  velocities  ahead  of  the  elastic 
o  S 

forerunner  wave  and  behind  the  forerunner  wave  respectively.  The  velocity  of  elastic 
forerunner  wave  itself  is  equal  to  ^  and  v^  is  the  particle  velocity  behind  the  shock 
wave.  From  initial  conditions,  v  =  0,  From  equations  (1.4)  and  (1.5)  the  shock  wave 
velocity  W^is  given  by  the  following  equation. 


.  Y  =  (&.f)Y 
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Usually,  vg  is  very  small  and 

Then 


This  is  the  velocity  of  the  loading  shock  wave  in  the  locking  material.  For 
example,  for  aluminum  locking  material  of  initial  density  pQ  =  1.5  gms/c.c., 
=  2.7  gms/c.c.,  Y  =  37,500  psi  (258.6  mpa)  and  aQ  =  125,000  psi  (862  mpa)  the  velocity 
of  the  loading  shock  wave  is  3354  feet  per  second  (1022  m/sec)  which  is  approximately 
20%  of  the  longitudinal  elastic  wave  velocity  of  the  solid  aluminium.  These  examples 
qualitatively  illustrate  the  slow  loading  velocities  in  a  locking  material  of  low  initial 
density.  However,  it  is  to  be  noted  that  the  shock  wave  velocity  in  a  solid  is  slower  than 
the  corresponding  longitudinal  elastic  wave  velocity.  Even  then,  the  shock  wave  velocity 
can  be  shown  to  be  slower  in  locking  materials  than  in  the  corresponding  solid. 


1.3  Experimental  Background 
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In  the  past,  experimental  investigations  have  been  conducted  on  several  locking 
materials.  In  these  investigations,  metals,  plastics,  graphite  and  ceramics  have  been 
used  as  parent  materials  to  produce  locking  materials.  Experiments  have  been  conducted 
by  using  light  gas  guns  and  flyer  plates.  In  these  experiments,  attempts  have  been  made 
to  simulate  conditions  of  one  dimensional  strain.  From  the  point  of  view  of  the  present 
study,  a  significant  result  from  these  studies  concerns  a  comparison  of  the  shock  wave 
velocities  in  locking  materials  with  the  corresponding  wave  velocities  in  solid  materials. 

In  most  of  these  studies,  a  two  wave  pattern  has  been  observed  in  locking  materials 
and  solids.  In  locking  materials,  a  forerunner  wave  carries  stresses  in  the  material 
before  the  phase  transition  from  the  initial  density  p  Q  to  the  locked  density  p^.  The 
second  slow  wave  corresponds  to  the  shock  wave  that  was  discussed  in  the  previous 
section.  This  shock  wave  carries  stresses  that  exceed  the  stresses  necessary  for  phase 
transition.  Solid  materials  also  exhibit  a  two  wave  pattern.  In  this  case,  the  stresses 
carried  by  the  forerunner  wave  correspond  to  the  yield  limit  and  the  elastic  behavior  of 
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the  solid  before  yield.  This  forerunner  wave  travels  at  a  velocity  corresponding  to  the 
longitudinal  elastic  wave  veiocity  in  the  isotropic  solid.  The  second  wave  carries  stresses 
that  exceed  the  elastic  limit.  The  experiments  confirm  the  fact  that  this  second  wave  in 
the  solid  is  much  faster  than  the  second  shock  wave  in  the  locking  material.  For 
example,  in  aluminum  locking  materials  of  initial  density  pQ  s  2.1  gms/c.c.,  the  second 
shock  wave  velocity  varied  from  0.7  to  1.29  mm/psec.  The  second  wave  velocity,  in  the 
corresponding  solid  was  in  the  range  of  4.62  mm  /y  sec.  Similarly,  the  first  wave 
velocities  in  the  locking  materials  were  smaller  than  the  first  wave  velocities  in  the 
corresponding  solid.  The  first  wave  velocities,  in  aluminum  locking  materials  of  initial 
density  pQ  =  2.1  gms/c.c.,  varied  in  the  range  1.6  mm/p  sec  to  2.0  mm/psec.  The 
corresponding  first  wave  in  the  solid  travelled  with  a  velocity  of  6.11  mm/psec. 

The  theoretical  foundation  and  the  results  of  the  experimental  studies  confirm  the 
potential  benefits  that  can  be  derived  from  locking  materials  when  they  are  used  as 
protective  structures. 

2.  RESEARCH  PLAN 

2.1  Problems  Associated  With  Engineering  Design 

In  the  point  of  view  of  engineering  designs,  the  results  of  this  project  should 
provide  tools  to  design  a  countermeasure  or  a  protective  structure  for  the  mine  blasts 
that  will  be  specified.  In  this  report,  however,  the  protective  structure  for  mine  blasts  is 
assumed  to  be  made  of  a  locking  material.  Specifically,  such  a  design  involves  the 
selection  of  a  parent  material  from  which  the  locking  material  is  produced,  the  thickness 
of  the  protective  structure,  other  geometrical  parameters  and  methods  of  joining  the 
countermeasure  to  the  structure.  In  addition,  the  selection  of  an  initial  density  for  the 
locking  material  and  the  specification  of  the  desired  microstructure  constitute  important 
aspects  of  the  design.  The  selected  initial  density,  the  specified  production  technique 
and  the  resulting  microstructure,  usually  determine  the  range  of  pressures  at  which  the 
phase  transition  takes  place  from  the  initial  density  to  the  locked  density.  The 
microstructure  is  also  responsible  for  the  elastic-plastic  stress-strain  behavior  and  the 
associated  yield  conditions. 

An  essen  'al  design  tool  for  the  selection  of  the  various  design  parameters,  is  a 
comp*  *er  sim»  -tion  that  is  capable  of  providing  the  transient  dynamic  response  of  the 
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locking  material  or  the  locking  material  and  structure  combinations.  The  specified  mine 
blasts  and  the  various  design  parameters  can  be  supplied  as  an  input  to  the  resulting 
computer  program.  The  output  can  be  obtained  as  deformations  and  stresses  at  various 
locations  as  a  function  of  time.  The  output  data  can  then  be  used 

(a)  to  check  if  the  selected  design  parameters  provide  the  desired  margin  of 
safety  for  various  failure  modes, 

(b)  to  iterate  or  to  modify  the  design  parameters  and 

(c)  to  establish  the  damage  tolerance  and  reliability  of  the  design. 

The  primary  objective  of  the  reported  investigations  is  to  develop  such  computer 

simulations  and  the  resulting  computer  programs.  To  date,  most  of  the  reported 

investigations  have  been  restricted  to  one  space  dimension  and  time.  The  results  have 

been  mainly  used  to  simulate  the  experimental  results.  The  programs  have  not  been  used 

to  study  the  effectiveness  of  the  locking  material  as  a  countermeasure  to  resist  mine 

blasts.  A  study  of  the  transient  dynamic  response  of  a  locking  material  that  involved  two 

u 

space  dimensions  and  time  has  been  reported  by  the  author  .  This  study  described  the 
impact  of  compactible  plates  under  axisymmetric  conditions.  The  constitutive 
relationship  that  was  used  for  the  compactible  plate  is  very  specialized  and  cannot  be 
easily  generalized  for  the  locking  materials  of  concern  here.  The  one  dimensional 
analysis  and  the  compactible  plate  analysis  have  used  the  finite  difference  technique 
developed  by  Von  Neuman, ®  and  Wilkins^.  Certain  difficulties  were  encountered  in 
the  application  of  the  Von  Neuman's  approach  to  locking  materials. 

(1)  Very  small  time  steps  were  needed  to  account  for  the  phase  transition  from 
the  initial  density  to  the  locked  density. 

(2)  An  accurate  computer  simulation  of  the  precursor  wave  was  very  difficult. 

(3)  Computational  efficiency  and  the  level  of  accuracy  decreased  when  two  or 
more  space  dimensions  were  considered.  Similar  problems  were  encountered 
in  considering  multiple  reflections. 

Other  investigators^  also  experienced  similar  difficulties  in  applying  the  Von 
Neumann's  technique  to  the  study  of  the  dynamic  response  of  locking  materials  and  phase 
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transition.  In  order  to  reduce  these  difficulties  and  to  develop  the  needed  computer 
simulation  techniques  for  transient  dynamic  response  of  locking  materials,  the  following 
research  plan  was  adapted. 


(1)  The  first  step  was  to  prepare  a  state-of-the-art  review  of  finite  difference 
methods  for  the  numerical  solution  of  hyperbolic  differential  equations. 

(2)  The  second  step  was  to  select  a  computationally  efficient  and  accurate  finite 
difference  method. 

(3)  The  selected  second  order  accurate  method  was  then  applied  to  study  the 
transient  response  of  locking  material  structures  under  the  conditions  of  one 
dimensional  strain.  The  purpose  of  the  analysis  was  to  demonstrate  the 
capabiltiy  of  locking  material  as  a  countermeasure  to  mine  blasts. 

(4)  Then,  a  computer  simulation  of  locking  material  was  developed  under 
conditions  of  axisymmetry  and  finite  deformations. 


3.  STATE-OF-THE-ART 


The  field  of  the  analysis  of  the  transient  dynamic  response  of  elastic-plastic¬ 
locking  materials  is  in  the  general  field  that  is  concerned  with  seeking  solutions  to 
nonlinear  hyperbolic  differential  equations.  For  both  linear  and  nonlinear  differentia1 
equations,  that  are  encountered  in  the  field  of  solid  mechanics,  very  few  analytical 
solutions  have  been  obtained.  As  a  consequence,  alternative  methods  of  solutions  have 
been  sought  for  many  practical  problems  of  solid  mechanics.  Integral  tranform 
techniques  have  been  very  effective  tools  in  solving  linear  elastodynamic  problems'*”'^. 
Primarily,  integral  transforms,  such  as  Laplace  tranforms,  have  been  used  to  remove 
time  as  an  independent  variable.  This  effectively  reduces  the  equation  to  an  elliptic 
type.  After  the  solutions  have  been  obtained  to  the  reduced  equation,  in  tranformed 
variables,  the  next  operation  is  to  obtain  the  inverse  tranform  to  return  the  solution  to 
the  time  domain.  One  of  the  particularly  successful  techniques  of  obtaining  the  inverse 

g 

transforms  for  elastodynamic  problems,  is  the  Cagniard-deHoop  technique  .  Basically, 

this  method  involves  the  modification  of  the  contour  of  integrations  that  are  encountered 

in  the  process  of  obtaining  inverse  transforms,  in  such  a  way  that  the  integrals  are 

rendered  to  a  form  with  known  solutions.  This  technique  has  been  successfully  applied  to 
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many  solid  mechanics  problems 

These  tranform  techniques  have  limitations  when  applied  to  nonlinear  problems  and 
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structures  with  finite  boundaries.  For  this  class  of  problems,  other  numerical  methods 

are  necessary.  Numerical  methods  that  use  the  method  of  characteristics,  have  been 

23 

primarily  restricted  to  the  class  of  one  space  dimensional  problems.  Butler  and 

2  L 

Clifton  have  reported  the  applications  to  two  space  dimensional  problems.  Most  of  the 
developments  in  the  last  two  decades  have  used  finite  difference  and  finite  element 
methods  to  solve  the  transient  dynamic  response  problems  in  solids. 

3.1  Finite  Difference  Methods 


The  research  activity  of  the  sixties  and  seventies  has  proved  that  the  finite 

element  methods  are  superior  to  other  numerical  methods  in  solving  elastostatic 

problems.  However,  the  same  superiority  does  not  apply  to  linear  or  nonlinear  transient 

dynamic  response  problems  in  solids.  In  fact,  the  finite  difference  methods  offer  a  viable 
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choice  to  the  solution  of  transient  dynamic  response  problems.  It  has  been  shown  that, 
for  selected  problems,  finite  difference  methods  provide  better  accuracy  and  computa¬ 
tional  efficiency  when  compared  with  the  solutions  obtained  by  using  finite  element 
methods.  The  rest  of  this  review  of  the  state-of-the-art  is  restricted  to  finite  difference 
methods  and  their  application  to  the  study  of  the  linear  and  nonlinear  transient  dynamic 
response  problems  in  solids. 

Early  work  in  the  field  of  the  application  of  the  finite  difference  methods  to 

transient  dynamic  response  of  solids,  consisted  of  the  application  of  a  finite  difference 

26 

scheme  developed  by  Von  Neumann  for  use  in  hydrodynamics.  The  Von  Neumann's 

scheme  is  an  explicit  finite  difference  approximation  of  the  time  derivatives  and  is  often 

27  28 

classified  as  "Leap  Frog"  scheme  .  A  classic  paper  in  this  field  is  by  Wilkins  to  solve 

29 

two  and  three-dimensional  transient  elastodynamic  problems.  Wilkins  has  also 

considered  ideal  plasticity  and  large  deformation  effects  in  his  studies. 

In  Wilkins  approach,  as  is  done  in  all  leap  frog  schemes,  the  dependent  variables 

are  staggered  in  space  and  time  to  satisfy  the  stability  requirements  and  maintain  second 

order  of  accuracy.  This  results  in  the  calculation  of  only  one  set  of  field  variables  at  any 

particular  mesh  point.  This  means  that  only  stresses  or  velocities  are  computed  at  a 

selected  mesh  point.  Furthermore,  such  a  staggered  scheme  sometimes  results  in 

computationally  induced  oscillations3®.  These  computationally  induced  oscillations  can 

be  minimized  by  using  an  artificial  viscosity  .  However,  the  use  of  the  artificial 

viscosity  results  in  a  finite  difference  scheme  that  is  not  optimally  stable3®.  The  time 

28  29 

steps  are  often  reduced  by  a  factor  of  three  in  two-dimensional  problems  ’  . 
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Variations  of  Wilkins'  approach  appear  throughout  the  literature  in  the  field  of 

24 

the  study  of  transient  elastodynamic  problems.  Clifton  has  proposed  a  different 

23 

approach  to  the  problem.  He  has  extended  the  difference  scheme  developed  by  Butler  , 

for  hydrodynamics,  to  study  two-dimensional  elastodynamic  problems.  This  finite 

difference  scheme  has  been  developed  by  formulating  integration  procedures  along 

bicharacteristics.  The  resulting  procedure  is  an  explicit  finite  difference  scheme. 

Unlike  the  leap  frog  scheme,  all  the  variables  are  calculated  at  all  mesh  points. 

Lax  and  Wendroff^  have  developed  a  second  order  accurate  scheme  for 

hydrodynamic  problems.  In  this  scheme,  all  the  dependent  variables  are  calculated  at  all 

mesh  points.  At  interior  points  Clifton's  scheme  and  Lax-Wendroff's  scheme  are 

identical.  Both  these  schemes,  however,  are  not  optimally  stable,  in  the  sense  of 
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Courant,  Friedrich  and  Levy,  for  more  than  one  space  dimension  .  This  leads  to 

computational  procedures  that  are  less  than  optimally  efficient.  Also  the  Lax-Wendroff 

schemes  involve  computation  of  squares  of  certain  matrices  and  result  in  a  complicated 
39 

algorithm.  Smith  has  applied  the  Lax-Wendroff  procedures  and  the  early  time  splitting 
40 

procedures  of  Strang  to  two-dimensional  problems.  Smith's  work  is  concerned  with  the 
comparison  of  the  relative  efficiency  of  the  schemes  and  is  not  a  complete  initial 
boundary  value  problem. 

However,  an  improved  version  of  Lax-Wendroff's  scheme,  that  has  the  potential  of 

providing  an  improved  computational  efficiency  and  accuracy,  can  be  applied  to  both 

linear  and  nonlinear  transient  dynamic  response  problems  in  solid  mechanics.  The 

improvement  can  be  incorporated  in  three  different  areas.  The  first  improvement  is  in 

40  40  41 

the  modification  of  the  scheme  to  make  it  optimally  stable  .  This  is  due  to  Strang  ' 

and  somewhat  similar  to  that  used  by  Smith.  The  second  improvement  is  to  improve  the 

computational  efficiency  by  providing  second  order  accuracy  at  every  other  time  step. 

38  42  42 

This  improvement  follows  the  developments  of  Gottlieb  ,  Gourlay  ,  Morris  and 

43  44 

Mitchell  .  The  third  improvement  is  in  the  incorporation  of  McCormack's  two  step 

procedure.  This  third  improvement  results  in  programming  simplicity  and  eliminates  the 

need  for  squaring  the  matrices  that  is  needed  in  Lax-Wendroff  or  Strang's  procedures. 

This  improvement  is  very  useful  in  solving  nonlinear  problems. 

An  improvement  in  the  implementation  of  the  finite  difference  schemes  is  the 

25 

concept  of  point  condition  codes  .  Attempts  at  improvement  of  osiilations,  overshoots 

45-51 

and  smearing  of  contact  discontinuities  are  also  possible 
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4.  ANALYSIS  UNDER  CONDITIONS  OF  ONE-DIMENSIONAL  STRAIN 
4.1  Development  of  the  Procedure 

The  proposed  second  order  accurate  method  of  analyzing  the  behavior  of  the 
locking  material  has  received  little  attention  when  applied  to  solid  mechanics  problems. 
The  efficiency  and  the  resulting  accuracy  of  the  method  have  been  studied  as  a  part  of 
another  sponsored  research  project  at  Georgia  Tech.  In  this  project,  one  and  two- 
dimensional  problems  of  transient  dynamic  analysis  have  been  investigated.  At  present, 
nonlinear  problems  have  been  studied.  In  the  previous  analyses,  the  nonlinear  behavior 
has  not  been  studied. 

In  this  section,  the  proposed  second  order  accurate  difference  scheme  has  been 
applied  to  the  transient  dynamic  response  of  locking  materials.  The  application  is 
restricted  to  one-dimensional  strain  in  the  x-direction.  The  numerical  results  are 
presented  for  the  specific  cases  of  loading  for  aluminum  locking  materials.  This  study  is 
undertaken  for  the  following  reasons. 

(1)  one-dimensional  problem  is  relatively  simple 

(2)  the  analysis  provides  an  understanding  of  how  the  proposed  second  order 
accurate  difference  schemes  can  handle  the  phase  transition 

(3)  the  effect  of  different  initial  densities  for  the  various  locking  materials  can 
be  qualitatively  understood  with  regard  to  a  one-dimensional  analysis. 

(4)  The  method  of  analysis  for  one-dimensional  problems  would  indicate  the 
possible  areas  of  difficulties  that  may  be  encountered  when  extended  to  two 
and  three-dimensional  problems  of  transient  dynamic  analysis  of  locking 
materials. 

The  one-dimensional  equations  of  motion  have  been  written  in  Lagrangian 
coordinates.  The  initial  positions  of  the  body  have  also  been  selected  as  the  Lagrangian 
coordinates.  Large  elastic-plastic  deformations  of  the  locking  material  have  been 
considered.  The  resulting  equations  are  as  follows: 

§4  (u)  =  A(p)  ||  {u}  (4.1) 


where 


{  u } T  =|  v(x,t),  p  (x,t),  a  x(x,t),  a  y(x,t)  j T 


(4.2) 
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In  these  equations,  it  has  been  assumed  that  the  stress  tensor  can  be  separated  into 
hydrostatic  pressure  and  stress  deviators.  The  hydrostatic  pressure  has  been  assumed  to 
be  related  to  the  changes  in  density  and  follow  a  locking  behavior  as  shown  in  Figure  1.3. 
The  stress  deviators  are  assumed  to  follow  an  elastic-ideal  plastic  behavior  with 
Von  Mises'  yield  condition. 

All  the  assumptions  of  the  preceding  paragraph  are  approximations.  However, 

these  approximations  have  been  used,  in  the  past,  to  express  the  mechanical  behavior  of 
3  u  52  52 

locking  materials  ’  ’  and  other  solids  at  very  high  stresses.  On  the  basis  of 

52 

experimental  results,  the  approximations  have  been  found  to  be  reasonable  .  In  the 
absence  of  any  available  constitutive  relationships  on  the  basis  of  second  Piola-Kirchhoff 
stresses  and  Green-Lagrange  strains,  this  paper  has  considered  similar  approximations. 
The  constitutive  relationships  have  been  written  in  terms  of  Cauchy  stresses,  particle 
velocity  gradients  and  density.  Then, 


a  =  -  p  ♦  S 

X  K  X 

(*.<0 

o  „  =  -  P  ♦  S 

(4.5) 

y  y 

p  =  fj(p) 

V 

o 

a. 

p  ■  f2<p) 

Pi  < 

P  <  Pt  i  ‘  (4*6) 

p  =  f3(p) 

»  J 
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The  quantities  f[t  f2  and  f3  are  as  foliows  (Figure  1.3): 


Similarly, 


W  =  P,  '  I  [(^)R  -  l] 

fj(p)  =  A(p/p£  -  1)  +  B(p/p^.  -  1)^ 


SxS-2Sy  =  .2S2  =  .iG£ 


with  the  yield  condition 


2  2  ?  2 
Sx  +  25y  -fY^o. 


<►.2  Numerical  Analysis 


A  numerical  integration  of  the  hyperbolic  partial  differential  equation  (4.1)  by  a 
procedure  similar  to  that  of  Lax  requires  that  the  field  variables  {  u  J  at  t  +  A  t  should  be 
calcualted  from  a  knowledge  of  the  field  variables  j  u  |  and  its  spatial  derivatives  at  t. 
To  a  second  order  accuracy, 


l“}t  +  ^  =  lijf*  +  jffj  At  *  +  0(^>3  •  C4.10) 

The  differential  equation  (4.1)  can  be  used  to  express  the  quantities  3u/3  t  and  3  ^u/3t^  in 
terms  of  the  spatial  derivatives  of  j  u  J.  Then, 

i-i*  •  ^ 

*  4"  *  [ff] { S } *  4r  •  «*•“> 
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The  equation  (4.11)  is  very  complicated.  It  involves  (A(p)l  and  multiples  of  first  and 
second  partial  derivatives  of  |  u  j  with  respect  to  x.  This  equation  can  be  considerably 
simplified  by  using  a  two  step  approximation  due  to  McCormack,  Richtmyer  or  Gottlieb. 
For  example,  McCormack  two  step  formulation  is  as  follows. 

,t 


f  ,  d  U  ) 

l“)  *  til  *  at  j  4' 


M'*4' =  Ml  ?!  ♦!“:!)  ♦*!£) 


(4.12) 


It  is  easy  to  verify  the  accuracy  of  (4.12)  by  expansion.  Now,  the  time  derivatives  on  the 
right  hand  sides  of  (4.12)  are  replaced  by  spatial  derivatives.  Thus, 

l“f  *  |;f  *  IA(»)I  a-f  *  At. 

#  # 

It  is  to  be  noted  that  (A(p)]  and  3y/dx  are  to  be  evaluated  by  using  the  star  values 
of  the  field  variables.  Also,  it  can  be  seen  that  the  formulation  described  in  the 
equations  (4.13)  is  much  simpler  than  the  formulation  described  in  (4.11).  To  complete 
the  finite  difference  formulation,  the  spatial  derivatives  of  j  u  J  and  j  u  J  are  expressed 
in  terms  of  their  spatial  finite  difference  approximations.  A  central,  forward  or 
backward  difference  is  used  depending  on  whether  the  point  under  consideration  is  an 
interior  point,  left  boundary  point  or  right  boundary  point. 

4.3  Stability  Requirements 

The  equations  (4.13)  can  be  written  in  the  form  of  finite  difference  operators. 

Then, 

L*  =rn+  k[AlAx 

Lj  =  ^IH  kt A 1  Ax 


m 


(4.14) 
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In  these  equations,  Ax  represents  the  spatial  finite  difference  approximations.  The 
quantity  I  is  the  identity  operator  and  At  has  been  replaced  by  k.  Then, 

|u|t+  At  =  <J-  +  LjL*)  |  u  j*  .  (4.15) 

The  operators  I,  Lj  and  L  ,  operating  on  J  u  j*  in  the  manner  shown  in  equation  (4.15) 
change  J  u  Jt  to  j  u  j*  +  This  is  an  explicit  scheme  and  is  subject  to  the  usual 
restrictions  on  k  =  At  to  maintain  the  stability  of  the  computational  scheme.  In 
particular,  it  can  be  shown  that  At  <Ax/c  where  c  is  the  fastest  of  the  local  elastic  wave 
velocities. 


4.4  Plasticity  Effects  and  Yield  Conditions 


The  numerical  method  that  has  been  used  in  this  paper  is  an  explicit  method.  The 

field  variables  at  t  +  At  depend  only  on  die  field  variables,  at  time  t.  The  particle 

velocities  and  strains  at  t  +  At  can  be  computed  from  the  current  field  variables  at  t. 

The  stresses  computed  at  t  +  At  may  violate  the  yield  condition  (4.9).  However,  these 

stresses  can  be  adjusted  along  appropriate  normals  to  the  yield  surface  back  to  the  yield 
52 

surface  .  The  already  calculated  velocities  and  strains  that  depend  only  on  the  field 
variables  at  t  remain  unchanged.  The  only  regions  where  iterative  calculations  are 
needed  are  at  the  boundaries. 


4.5  Discussion  of  Results 

In  order  to  achieve  the  same  accuracy  as  a  first  order  accurate  method,  a  time  step 
equal  to  the  square  root  of  A  t  that  is  used  in  the  first  order  method  is  needed.  The 
quantity  At  is  usually  less  than  one  and  hence  the  second  order  accurate  method  uses  a 
larger  time  step.  Hence,  a  larger  At  can  be  chosen  to  obtain  the  same  order  of  accuracy 
as  a  first  order  method.  This  leads  to  a  smaller  number  of  the  finite  difference  cells  and 
increased  computational  efficiency.  The  developed  computer  program  has  been  used  to 
examine  the  efficiency  of  locking  materials  in  attenuating  stresses  and  protecting 
structures.  In  this  analysis,  specified  items  include  the  peak  impact  stress  and  the 
unloading  pattern  of  the  stress  wave  impinging  on  the  locking  materials.  The  locking 
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density  is  usually  fixed  for  a  given  locking  material.  The  results  of  the  numerical 
analysis  can  be  used  to  explain  the  peak  stress  and  the  stress  distribution  for  various 
times  after  impact  and  unloading  at  the  boundary. 

In  the  numerical  example,  aluminum  locking  materials  have  been  selected.  Some  of 
the  properties  of  aluminum  locking  materials  have  been  extensively  investigated3.  A 
pressure  density  relationship  as  shown  in  Figure  1.3  has  been  assumed.  The  locking 
density  is  equal  to  the  2.72  gms/c.c.  Various  initial  densities  have  been  considered. 
Results  have  been  presented  for  pQ  *  1.39  and  2.1  gms/c.c  here.  There  are  three  distinct 
branches  of  the  pressure-density  relationship.  Equations  for  these  branches  are  the  same 
as  in  equations  (4.6)  and  (4.7). 

The  quantity  is  the  intersection  of  the  p  -p  curve  with  p  axis  as  shown  in  Figure 
1.3.  The  appropriate  constants  are  selected  from  reference  3.  For  all  P>  Pg  the 
unloading  is  assumed  to  follow  the  slope  of  the  solid  p  -  p  curve  as  shown  in  Figure  1.3. 

Two  types  of  loading  have  been  considered.  The  first  type  of  loading  consists  of  a 
step  loading  followed  by  step  unloading.  It  has  been  assumed  that  a  peak  compressive 
stress  of  oQ  is  applied  at  time  t  =  0  to  the  left  boundary  of  the  slab  of  a  locking  material. 
For  purposes  of  illustration,  a  value  of  -  aQ  =  100,000  psi  (689.5  MPa)  has  been  assumed. 
The  applied  stress  is  reduced  to  zero  by  step  unloading  at  time  t  *  0.2  ysec.  This  loading 
pattern  is  illustrated  in  Figure  4.1  and  will  be  called  loading  pattern  'a'.  The  second  type 
of  loading  consists  of  a  step  loading  of  a  compressive  stress  of  magnitude  -oQ  at  t  =  0 
followed  by  an  exponential  unloading  at  the  left  boundary.  For  purposes  of  illustration, 
-Oq  has  been  assumed  to  be  equal  to  100,000  psi  (689.5  MPa).  The  applied  stress  is 
assumed  to  be  maintained  at  a  ~  -Og  for  a  duration  of  0.02  y  sec.  Then  the  exponential 
decay  of  the  applied  load  at  the  left  boundary  decreases  the  magnitude  of  Oq  to  .05  <Jq  at 
0.2  usees.  This  unloading  pattern  is  illustrated  in  Figure  4.2  and  will  be  called  loading  V. 

First,  a  solid  material  slab  of  initial  thickness  0.5  inches  (12.7  mm)  has  been 
considered.  The  thickness  has  been  divided  into  100  cells.  In  this  case  the  solid  density 
is  =  2.72  gms/c.c.  The  transient  response  of  this  slab  to  the  loading  pattern  'a'  has 
been  studied  by  using  the  developed  computer  program  and  the  results  are  illustrated  in 
Figure  4.3.  This  figure  is  a  plot  of  the  stress  a  x  as  a  function  of  the  thickness  at  two 
different  instants  of  time  t  =  0.2 1  y  secs  and  t  =  0.68  y  secs.  As  seen  in  the  figure,  the 
stress  at  a  distance  of  0.14  inch  (3.6  mm)  from  the  initial  left  boundary  is  still 
approximately  100,000  psi  (689.5  MPa).  This  is  equal  to  the  applied  peak  stress.  No 
significant  attenuation  has  taken  place  during  the  travel  of  the  stress  wave  through  the 
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thickness  equal  to  0.14  inch  (3.6  mm). 

Next,  a  locking  material  of  initial  density  pQ  =  2.1  gms/c.c.  and  locking  density 
=  2.72  gms.  c.c.  has  been  considered.  The  thickness  and  division  into  cells  are 
Identical  to  those  for  the  solid.  The  transient  response  to  the  loading  pattern  'a'  has  been 
computed  and  illustrated  in  Figure  4.4.  The  stress  distribution  as  a  function  of  the 
distance  from  the  left  boundary  has  been  illustrated  for  t  =  0.13  p  sec,  0.54  y  sec, 
1.15  y  sec,  and  2.22  y  sec.  It  can  be  seen  that  the  peak  stress  has  reduced  by 
approximately  35%  over  a  distance  of  0.09  inch  (2.3  mm)  from  the  left  boundary.  In  this 
figure,  the  elastic  forerunner  wave  can  also  be  seen.  Similar  stress  distributions  are 
observed  for  locking  materials  of  initial  density  p  =  1.818  gms/c.c.,  1.604  gms/c.c.  and 
1.39  gms/c.c.  If  these  locking  materials  are  considered  as  countermeasures  in  front  of  a 
given  structure,  the  peak  impact  stress  has  been  reduced  by  an  amount  as  much  as  50% 
when  the  shock  wave  has  propagated  a  distance  of  0.013  inch  (0.33  mm)  (Figure  4.5).  It 
can  also  be  observed  that  the  locking  material  of  lower  initial  density  has  the  potential 
of  attenuating  the  impact  peak  stress  by  a  larger  percentage  when  compared  with  the 
locking  material  of  higher  initial  density. 

As  a  next  step,  the  loading  pattern  V  has  been  considered.  The  results  of  the  study 
of  transient  response  through  a  solid  has  been  illustrated  in  Figure  4.6.  The  results  are 
similar  to  those  for  the  loading  pattern  'a*.  The  peak  impact  stress  has  not  attenuated 
during  the  passage  of  the  stress  wave  over  a  distance  of  0.15  inch  (3.8  mm).  Similar 
stress  distribution  for  loading  pattern  V  has  been  illustrated  in  Figures  4.7  and  4.8  for 
different  locking  materials  of  initial  density  pQ  =  2.1  gms/c.c.  and  1.39  gms/c.c.  The 
Figure  4.7  is  for  a  locking  material  with  initial  density  2.1  gms/c.c.  Over  a  shock  wave 
traverse  of  .045  inch  the  peak  stress  has  been  reduced  by  60%.  However,  the  locking 
material  with  initial  density  of  1.39  gms/c.c.  can  attenuate  the  peak  stress  by  almost 
90%,  during  the  shock  traverse  of  0.05  inch  (1.3  mm).  This  is  illustrated  in  Figure  4.8. 
The  unloading  waves  and  the  elastic  forerunners  can  be  identified  in  all  these  figures.  In 
Figure  4.9,  the  effect  of  the  plasticity  on  the  transient  dynamic  response  is  shown.  A 
yield  stress  of  80,000  psi  is  assumed  for  purposes  of  illustration. 

5.  ELASTIC  -  PLASTIC  -  LOCKING  MATERIALS  UNDER  CONDITIONS  OF  AXISYMMETRV 
5.1  Introduction 

In  this  section,  the  numerical  scheme  that  has  been  discussed  in  Sections  3  and  4 
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has  been  modified  and  applied  to  study  the  transient  dynamic  response  of  elastic-plastic 

locking  materials  under  conditions  of  axisymmetry.  A  cylindrical  polar  coordinate 

system  r,  0,z  has  been  used.  Initial  positions  of  various  material  points  have  been  used  as 

Lagrangian  coordinates.  The  problem  has  been  formulated  as  a  finite  deformation 

problem  with  Cauchy  stresses  that  are  defined  in  a  deformed  coordinate  system.  In  order 

to  represent  the  equations  in  a  form  that  is  suitable  for  the  use  of  Gottlieb 
uu  40  41 

MacCormack  -Strang  ’  type  of  scheme,  the  constitutive  relationships  have  been 
used  in  a  rate  form.  The  Jaumann  stress^  rate  (or  the  corotational  stress  rate)  has  been 
used  to  satisfy  the  principle  of  objectivity. 

5.2  Governing  Equations 

The  governing  equations  are  then  written  as  follows. 

Kinematic  equations: 


Equations  of  motion : 


i  =  w 

(5.1) 

t  =  u 

(5.2) 

—1  a 

ti 

•> 

3%z  1 

■  ■  ^  «• 

3Tzr  1  Tzr 

(5.3) 

a  z  p 

3r  p  r 

1 

11=  — ■ 

aV  i 

1  X  — 

3orr  1  arr-°06 

(5.4) 

U  P 

az  p 

a  r  +  p  r 

Constitutive  equations: 


*»  -  *  J®  %  *  1®  If  *  <4W®  <”> 


v  -  <p£-!®B*  *  J«l?»  *$-!«?♦  <.«-!?>  <»> 
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Jee 


(5.7) 


zr 


-  3u  -  3w  1  / 
Qfi  ♦  Gr.  ♦  7(0, 


3  r  T  2  zz  rr 


(5.8) 


Continuity  equation: 


P 


(5.9) 


It  has  been  assumed  that 


p  =  p(p)  (5.10) 

follows  the  equations  similar  to  (4.6)  that  represent  a  locking  behavior.  A  Von  Mises 
yield  condition  and  ideal  plasticity  have  been  assumed. 


*azz  +  p)2  +  *arr  +  P*2  +  ^CT00  +  p^2  +  2lzr2-l  ^ 


The  notation  used  in  the  equations  are  as  follows: 


z,r 


w 

u 


P 

P 

G 

Y 

L 

R 


current  positions  of  the  coordinate 
velocity  in  z-direction 
velocity  in  r-direction 
normal  stress  in  z-direction 

normal  stress  in  r-direction 

normal  stress  in  -direction 

shear  stress 

density 

hydrostatic  pressure 
modulus  of  rigidity 
yield  stress  in  simple  tension 
length,  in  the  z-direction 
radius  of  the  structure 


'1 
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"*  ',**u°ns  ,iU  ,0  M>  <*"  <*  "writtw  in  the  following  form 

!“)=  ju},r.Cc](u) 

where 


(5.12) 


r,  w 


’u’az2’V 


aee>T2f»pjT  • 


The  matrices  and  [c]  are  not  constants 

elastic-piastic-iocking  materials.  Specifically 


in  the  finite  deformation  problem  of 


0 

0 

0 


0 

0 

0 


dE  +  *G 
dp  +  3 

dE  _  2G 
dp  3 

dE  2G 
dp  ”  3 


-t 


zr 


zr 


0  G+={a  -  o  ) 

2  zz  rr 


-P 


0  0  0  0  0 

0  0  0  0  0 

i/p  o  o  o  o 

0  0  0  I/p  0 

0  0  0  0  0 

ooooo 
0  0  0  0  0 

ooooo 

o  0  0  Q  0 


(5.13) 
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0 
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0 

0 

0 
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0 

Q 

0 

0 

0 

0 

0 

0 

0 

Tzr 

odE  2C 
P  dp  -  T 

0 

0 

-T 

0dE  + 

zr 

Pdp  *  T 

0 

0 

0 

G-|(o 

2  zz 

pg£  - 

p  dp  T 

0 

0 

-0)0 

rr 

0 

0 

0 

-  p 

o  0  0  o  "o 

OOOOO 
0  0  0  1/P  o 

0  i/p  0  o  o 

0  o  0  o  o 

0  o  o  o  o 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 


(5.14) 
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1  0 
0  1 
0  0 


0  0 

0  <P^-f>/ r 

0  (p*-f)/r 

0  (pSE  .  f  )/r 
0  0 

0  -p/r 


0  0  0  0  0 

0  0  0  0  0 

0  0  0  1/pr  0 

0  1/pr  -1/pr  0  Q 
0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 


(5.15) 


5.3  Numerical  Method 


For  the  hyperbolic  differential  equation  (5.12),  the  Lax-Wendroff  second  order 
operator  can  be  written  as  follows. 


(ul  '  *  4*  .  Lzr{u!* 


(5.16) 


where 


at2 

* 


r  z 


Lzr  =  [l]  +  it(CA]A2  +  CB]\  +  [C]) 

[A]2  4ZZ  *  [®]2  arr  *  <M  [B]  *  [BIA]> 

k. 

[A][A}U^Z  *  [B]  [  A)'u  4r  Ar 
[A]’U[AK  *  W.W  4r  *  [A]’„[c]  4o  4Z 


.  ([A][C].  [C]  [A])  4; 
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time  step  is  restricted  by  the  following  equation. 


« 


At 

min(Az.Ar) 


max(X  a>Xd)  <  — 

A  B  Vs 


where  X^,  Xg  are  the  eigenvalues  of  A  and  B  respectively.  In  an  optimally  stable 

scheme,  the  right  hand  side,  which  is  also  often  called  Courant,  Friedrichs  and  Levy 

number  (CFL),  should  be  1,  instead  of  1/V8*  Furthermore  the  operator  (5.17)  is 

complicated  and  the  increased  number  of  computations,  in  many  cases,  negates  any 

computational  efficiency  obtained  by  resorting  to  second  order  accurate  methods.  The 

40  41 

Lax-Wendroff  operator  can  be  replaced  by  a  modified  Strang's  ’  operator  to  make  the 
Lax-Wendroff's  scheme  optimally  stable.  The  modification  of  Strang's  scheme  involves 
replacing  the  central  difference  operators  of  Strang  by  contour  differences  and 
accomodating  the  axisymmetric  equations.  Then, 


{u}t+At  =  LAt^2  LAt/2  LAt/2  L  At^2{  u}1  (5.19) 

where 

L“/2  =  [']*  t  <0K  *  [°k> 

*  <!>][>]*  [°][A] 

*  *  [A]’u[A]4r  *  [AM°]  O1  4Z 

*  <[D]2  *  [A][D]’U4Z  *  [D}u[Ak  *  [D]'U[D]vj  (5-20) 
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Lr‘,/2  =  [']*¥<[*  K  *M> 

*  ^{["Pv  *  <[b][e]  *  [>][>]♦  [b]H'u  \ 

*[B}u[B]4r  *  [B]-u[EK)4r  *  <[E][EMB][E]vA 


[E]’IB]4r  *  [E1U[EK> 


(5.21) 


[D].[E].[c] 


(5.22) 


Even  though  the  Strang's  type  of  operator  (5.20)  is  optimally  stable,  the 

computational  efficiency  can  be  substantially  increased  by  restoring  the  second  order 

accuracy  at  every  other  time  step  by  following  the  procedures  similar  to  those  of 
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Courlay,  Morris  and  Mitchell  ’  .  Then 


L  t  +  2it  ,  L»t  .At  l4.  tAt{  ,  . 

rz  r  z  z  r  1  ‘ 


(5.23) 


The  finite  difference  operators  (5.23),  (5.20)  and  (5.21)  still  contain  multiplication  of 

matrices,  such  as,  [  A  }  (B),  |  A) ,  1  A)  etc.  The  operation  can  be  simplified  by  using 

44  u 

two  step  procedures  .  Thus,  lor  example 


L  *  =  (I]+At(AjA  +AtlDl 
2  2 

Lz  =  i(fll  +  Lz*)  ♦  (lA]#Az  +  L  D  ]*) 


(5.24) 
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*  #  # 
where  I  A 1  and  [  D  ]  are  evaluated  by  using  the  value  of  {  u }  .  A  similar  expression 

can  be  written  for  Lr  in  r-direction.  Suitable  spatial  contour  difference  forms  that 

maintain  second  order  accuracy  and  are  suitable  for  use  in  equations  of  the  type  (5.24) 

* 

have  been  developed  at  Georgia  Institute  of  Technology  .  These  are  different  from  those 
used  in  references  (28,  29) 

For  an  interior  point,  a  typical  member  V  of  the  field  vector  j  u  j  is  evaluated  as 
follows. 


V(1)  =  U*  +  At(BAr  ♦  E)Ul 

V<2)  *  5(U*  *  V  *  *  E(l)»v(l) 


V(3)  =  V(2)  +  At(A(2)Az  +  D(2))V(2) 


V(4)  =  I(V( 2)  +  V(3)}  +  (A(3)7z  +  °(3))V(3) 


V(5)  =  V(4)  +  At(A(4)  Az  +  D(4)}  V(4) 


(5.25) 


V(6)  =  I(V(4)  +  W  +  2At(A(5)7z  *  °(5))V(5) 


V(7)  =  V(6)  +  A  t(B(6)  Ar  +  E(6))  V(6) 


V(8)  =  2  (V(6)  +  W(7))  +  IAt(B(7)7r  +  E(7))W{7) 


yt  +  2  At  =  y 


(8) 


In  this  equation  Af  and  7r  represent  the  contour  difference  operators  that  represent  the 
equivalent  forward  and  backward  difference  operators,  that  are  needed  in  the  two  step 


#Chen,  H.  P.,  Ph.D.  Thesis 


25 


procedures.  At  the  boundaries,  these  equations  have  been  modified  to  include  the 
appropriate  boundary  conditions. 

5.4  Numerical  Results  and  Discussion 


The  material  considered  for  the  detailed  numerical  analysis  has  the  following 
properties: 


Solid  density,  p^  =  2.72  gms/cc 
Young's  Modulus,  E  =  1.1  x  10^  psi  (75845  MPa) 

Poisson's  ratio,  v  =0.3 

In  the  following  analysis,  three  different  types  of  material  behavior  have  been 
considered.  These  are:  (i)  elastic  behavior,  (ii)  elastic-plastic  behavior  and  (iii)  elastic- 
plastic  behavior  in  a  locking  material.  For  the  study,  a  yield  stress  Y  =  20,000  psi 
(137.9  MPa)  has  been  assumed.  For  the  locking  material,  the  initial  density  is  chosen  to 
be  2.1  gms/cc.  There  are  three  distinct  branches  of  the  pressure-density  relationship. 
These  are  shown  in  Figure  5.1.  For  this  analysis,  first  two  branches  are  assumed  to  be 
straight  line  segments.  The  values  chosen  for  pQ,  pj,  Pj  etc.  (figure  5.1),  for  the  present 
numerical  studies,  are 

PQ  =  2.1  gms/cc  Pj  =  6600  psi  (45.5  MPa) 

P|  =  2.11  gms/cc  pg  =  66600  psi  (459.2  MPa) 


and 


P4  =  2.72  gms/cc  • 

The  equations  of  the  lines  corresponding  to  the  three  branches  of  the  pressure-density 
relations  are 


(  p  -  P0) 


P  <  Pj 


P 


Pl-Po 
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PQ-  Pi 

P  -  Pi  *  5-^7  to  -  Pi>  P,  <  P  <  0t 
o  1 

p  =  A  (p/PG  -  1)  +  B  (p/pg  -  l)2  p  >  P4  . 

The  constants  A,  B,  pG  are  selected  from  reference  3. 

Two  different  time  variations  of  loading  are  considered.  The  first  type  is  a  step 

loading  (without  unloading)  as  shown  in  Figure  3.2.  The  second  type  is  also  a  step  loading 

but  with  an  unloading  at  a  time  t  =  t.  This  is  shown  in  Figure  4.1.  The  magnitude  of  the 

peak  compressive  stress  0  has  been  assumed  to  be  100,000  psi  (689.5  MPa).  The 

geometry  of  the  structure  and  the  corresponding  applied  loading  are  shown  in  Figures 

5.3a  and  5.3b.  The  boundary  conditions  at  the  unloaded  edges  are  also  shown  in  the  same 

figures.  These  boundary  conditions  are  the  same  for  both  the  cases  of  applied  loading. 

Physically,  these  boundaries  can  be  considered  to  be  "frictionless-rigid".  This  means 

that,  on  these  edges,  the  velocities  in  the  direction  normal  to  the  plane  of  edges  and  the 

shear  stresses  are  zero.  Referring  to  figures  5.3a  and  5.3b,  these  boundary  conditions  are 

expressed  as  w  =  0,  r  =  0  at  the  boundary  which  can  be  denoted  by  z(t  =  0)  =  L  and 
zr 

u  =  0,  t  =  0  at  the  boundary  denoted  by  Kt  =  0)  =  R.  On  the  left  boundary,  represented 
zr 

by  zft  =  0)  3  0,  there  are  two  types  of  applied  loading.  One  type  corresponds  to  the  load 
on  the  entire  edge  as  shown  in  Figure  5.3a  and  the  other  corresponds  to  a  partial  load  on 
the  edge  as  shown  in  Figure  5.3b.  There  are  two  cases  of  loading  type  shown  in  Figure 
5.3a.  These  cases  are  denoted  as  case  5a  and  case  5b.  In  both  these  cases,  the  loading, 
the  geometry  and  the  boundary  conditions  are  such  that  a  one-dimensional  strain 
condition  is  created  in  the  axisymmetric  problem  that  consists  of  two  space  dimensions 
and  time.  The  problems  under  these  cases  are  solved  in  order  to  simulate  a  one¬ 
dimensional  problem  of  transient  dynamic  analysis  that  can  be  compared  with  known 
solutions.  The  difference  between  case  5a  and  case  5b  is  that  case  5a  has  only  a  step 
loading  while  case  5b  has  a  step  loading  followed  by  a  step  unloading.  The  case  5b  can 
also  be  called  a  pulse  loading. 

Similar  to  the  above  cases,  there  are  also  two  different  cases  under  the  partial 
loading  condition  shown  in  Figure  5.3b.  These  cases  are  identified  as  case  5c  and  case 
5d.  The  partial  loading  creates  two-dimensional  conditions.  Again,  the  difference 
between  case  5c  and  case  5d  is  that  case  5c  represents  a  step  loading  and  case  5d 
represents  a  pulse  loading. 
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The  dimensions  pertaining  to  the  problem  are  as  follows  (Figure  5.3) 


L  =  3.2  in  (81  mm),  R  =  0.8  in  (20  mm),  Rj  =  0.32  in  (8  mm) 


In  order  to  apply  the  finite  difference  scheme,  the  region  bounded  Internally  by  the  edges 
z(t  =  0)  =  0,  L  and  r(t  =  0)  =  0,  R  is  divided  into  WO  cells.  The  dimensions  of  these  cells 
are:  Az  =  Ar  =  0.08  in  (2  mm).  This  mesh  pattern  is  used  for  all  the  cases.  Since  the 
finite  difference  scheme  proposed  in  this  section  is  optimally  stable,  the  CFL  number  is 
chosen  to  be  1.  However,  the  time  step  At  is  not  always  a  constant.  All  the  plots  display 
the  variation  of  the  stress  a  zz  with  repsect  to  the  deformed  z-coordinate.  Also,  the 
symbols  E,  E-P  and  E-P-L  in  the  figures  refer  to  the  elastic  solid,  elastic-plastic  solid 
and  elastic-plastic-locking  material  respectively. 

As  a  first  step,  some  typical  results  of  the  stress  variations  obtained  by  using  a  leap 
frog  scheme  (programed  at  Georgia  Institute  of  Technology)  with  linear  or  quadratic 
artificial  viscosities  are  compared  with  those  obtained  by  the  second  order  accurate 
method.  These  are  shown  in  Figures  5.4a,  5.4b  and  5.4c.  In  these  figures,  step  loading 
condition  of  case  5a  has  been  assumed  and  only  elastic  behavior  is  considered.  The 
addition  of  artificial  viscosity  in  a  leap  frog  scheme  results  in  a  scheme  that  is  no  longer 
optimally  stable.  Suitable  CFL  numbers  must  be  chosen  in  these  numerical  calculations. 
In  Figure  5.4a  the  descriptions  of  the  plots  correspond  to  the  CFL  numbers  and  times  t  as 
given  in  the  following  table. 

Method  used  CFL  number  Elapsed  time  t  after 

the  application  of  load 

present  method  1.0  9.25  ysecs 

leap  frog  scheme  0.3334  9.26  ysecs 

with  linear  viscosity 


leap  frog  scheme  0.3334 

with  quadratic  viscosity 


9.51  ysecs 


It  can  be  seen  from  Figure  5.4a  that  the  a  -  z  variation  by  the  present  second  order 
accurate  method  has  a  very  sharp  wave  front,  and  indicates  very  little  oscillation  behind 
the  shock  wave.  A  modest  amount  of  overshooting  of  the  shock  front  can  also  be 
observed.  The  curve  obtained  by  leap  frog  scheme  with  linear  viscosity  has  a  very  wide 
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transition  layer  instead  of  a  steep  shock  front,  although  it  has  no  overshooting  and 
oscillations  behind  the  shock  wave.  In  the  plot  that  represents  the  results  of  the  leap 
frog  scheme  with  quadratic  viscosity,  oscillations  behind  the  shock  wave  are  significant. 
The  solution  represents  a  problem  with  one-dimensional  strain  and  hence  no  physically 
significant  oscillations  are  present.  All  oscillations  are  computationally  induced.  A 
larger  overshooting  is  observed  when  compared  with  that  of  the  second  order  accurate 
method.  The  shock  wave  front  is  also  not  steep.  The  results  of  using  higher  CFL 
numbers  in  the  leap  frog  scheme  with  artificial  viscosity  are  shown  in  Figures  5.4b  and 
5.4c.  For  the  results  shown  in  Figure  5.4b,  the  CFL  number  used  is  0.66  for  both  the 
linear  and  quadratic  viscosity  cases.  Obviously,  the  results  indicate  instability. 
Oscillations  behind  the  wave  front,  large  overshooting  and  wide  transition  shock  layers 
can  be  observed.  In  Figure  5.4c,  the  values  of  CFL  =  0.7  and  0.9  are  used  in  the  leap  frog 
scheme  with  linear  and  quadratic  viscosities  respectively.  These  results  are  also  not 
satisfactory  because  of  the  numerical  instability.  From  this  elementary  analysis,  it  can 
be  concluded  that  the  second  order  accurate  finite  difference  method  suggested  in  this 
report  is  a  desirable  choice  for  using  larger  time  steps. 

Elastic,  elastic-plastic  and  elastic-plastic  locking  behavior  are  compared  in  Figures 

5.5  and  5.6.  In  both  these  figures,  the  stress  c  has  been  plotted  against  the  current  z 

zz 

coordiantes  that  were  originally  located  along  Kt  =  0)  =  0.  The  Figure  5.5  corresponds  to 
approximately  7  usees  after  loading  at  die  boundary  z(t  =  0)  =  0.  The  corresponding  time 
for  Figure  5.6  is  approxiamtely  equal  to  12  usees.  The  specific  times  for  the  elastic 
body,  elastic-plastic  and  elastic-plastic-locking  bodies  in  the  Figure  5.5  are  6.61  usees, 
6.60  usees  and  6.74  Usees.  Similar  times  in  Figure  5.6  are  11.09 usees  for  the  elastic 
body,  11.86  usees  for  the  elastic-plastic  body  and  11.96  usees  for  the  elastic-plastic 
locking  material.  The  calculated  values  were  available  at  these  different  times  because 
of  the  different  values  of  At.  The  three  different  wave  fronts  can  be  seen  for  the 
elastic-plastic  locking  material.  The  single  wave  front  for  die  elastic  body  and  a  two 
wave  front  system  for  elastic-plastic  body  can  also  be  seen.  There  are  no  oscillations 
behind  the  wave  front  because  of  the  simulation  of  one-dimensional  strain  conditions  by 
a  combination  of  loading,  geometry  and  boundary  conditions.  The  velocity  of  the  locking 
wave  is  the  slowest  of  the  waves. 

In  Figure  5.7,  die  results  for  die  case  of  a  pulse  loading  are  shown.  The  loading 
consists  of  a  step  loading  followed  by  a  step  unloading.  Again  the  loading,  geometry  and 
the  boundary  conditions  simulate  the  one-dimensional  strain  conditions.  The  applied  load 


29 


is  set  to  zero  at  t  =  t.  The  value  of  t  =  1.33  p  secs  for  the  elastic  body,  t  =  1.33  p  secs 
for  the  elastic-plastic  body  and  t  =1.42  p  secs  for  the  elastic-plastic  locking  material 
have  been  selected.  The  figure  clearly  indicates  that  the  attenuation  of  peak  stresses  is 
maximum  in  the  locking  material.  A  reduction  of  peak  stress  is  also  seen  in  the  elastic- 
plastic  body  without  locking.  The  reduction,  however,  is  not  significant  in  comparison  to 
the  locking  material.  In  the  elastic  body  the  peak  stresses  are  not  reduced.  These 
results  support  the  fact  that  locking  materials  are  effective  stress  attenuators. 

The  loading  for  cases  5c  and  5d  correspond  to  a  partial  loading  on  the  surface 
z(t  =  0)  =  0.  This  corresponds  to  a  two-dimensional  (r  -  z)  problem.  The  variation  of  ozz 
against  z  are  shown  in  Figures  5.8  to  5.10.  These  are  for  a  step  loading  without 
unloading.  The  two  Figures  5.8  and  5.9  correspond  to  lines  r(t  =  0)  =  0  and  0.56  in 
(14.2  mm).  It  is  observed  that  the  elastic-plastic  locking  material  displays  a  steeper 
wave  front  in  comparison  with  an  elastic-plastic  material.  The  oscillations  behind  the 
elastic  wave,  due  to  boundary  effects  in  this  two-dimensional  problem,  can  be  seen.  The 
stress  magnitudes  in  Figure  5.9  are  small.  This  is  due  to  the  fact  that  the  region  is  away 
from  the  area  of  the  load  application.  The  variations  of  a  zz  with  z,  for  r(t  =  0)  =  0  and 
the  locking  material,  are  shown  in  Figure  5.10  for  increasing  values  of  times.  The 
development  of  three  distinct  wave  fronts  can  be  seen. 

The  Figures  5.11  and  5.12  correspond  to  pulse  loadings  on  a  part  of  the  surface 
specified  by  z(t  =  0)  =  0.  Again,  two-dimensional  effects  are  present.  The  Figure  5.11  is 
for  r(t  =  0)  =  0.  The  Figure  5.12  is  for  r(t  =  0)  =  0.56  in  (14.2  mm)  which  is  located  away 
from  the  loading  region.  Again,  observation  regarding  the  attenuation  of  peak  stresses 
can  be  made. 

6.  CONCLUSIONS 

In  this  report,  it  has  been  shown  that  a  locking  material  shield  has  the  potential  of 
protecting  a  given  structure  from  dynamic  loads  that  have  a  time  history  similar  to  that 
of  mine  blasts.  This  conclusion  is  based  on  the  analysis  under  conditions  of  one¬ 
dimensional  strain  and  the  studies  under  conditions  of  axisymmetry.  In  order  to  provide  a 
tool  for  designing  the  locking  material  shield,  a  second  order  accurate,  computationally 
efficient  numerical  scheme  has  been  developed.  The  numerical  scheme  is  capable  of 
considering  large  deformations  and  locking  phase  transition  and  provide  accurate  results 
in  comparison  to  the  present  state-of-the-art.  Point  condition  codes  are  also  developed 
for  locking  material-structure  combinations.  The  accuracy  and  the  reliability  of  the 
resulting  computer  program  have  been  checked  by  comparing  the  results  with  known 
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solutions. 

To  follow  the  present  work,  some  additional  tasks  are  recommended.  One  such  task 
is  to  improve  the  accuracy  and  the  computational  efficiency  during  phase  transition  and 
unloading.  Furthermore,  the  shockwave  can  be  steepened  by  employing  new  techniques 
such  as  the  method  of  artificial  compression.  Another  area  of  suggested  work  is  the  non¬ 
ideal  plasticity.  A  very  important  additional  task  will  be  to  investigate  the  accuracy  of 
methods  such  as  the  state  space  approach  with  its  inherent  exactness  and  simplicity  of 
reducing  three-dimensional  problems  to  those  involving  only  two-dimensions. 
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Fig. 1.1  A  Simple  Locking  Material 


Fig. 1.2  An  Elastic  Locking  Solid 
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